#include "borwein.h"



void calculationPiSeqBorwein() {

	// valor da iteracao atual
	int k = 0;

	// variaveis gmp
	mpf_t a, y, t1, t2, t3, t4, pi_anterior, pi_atual;

	// inicializacao de a, y e temporarios pelo gmp
	mpf_init2(a, 10000000);
	mpf_init2(y, 10000000);
	mpf_init2(t1, 10000000);
	mpf_init2(t2, 10000000);
	mpf_init2(t3, 10000000);
	mpf_init2(t4, 10000000);
	mpf_init2(pi_anterior, 10000000);
	mpf_init2(pi_atual, 10000000);


	// valores iniciais
	// a[0] = 6 - 4 * raiz(2)
	mpf_sqrt_ui(a, 2);		// a[0] = raiz(2)
	mpf_mul_ui(a, a, 4);	// a[0] = raiz(2) * 4
	mpf_ui_sub(a, 6, a);	// a[0] = 6 - a

	// y[0] = raiz(2) - 1
	mpf_sqrt_ui(y, 2);		// y[0] = raiz(2)
	mpf_sub_ui(y, y, 1);	// y[0] = y[0] - 1


	// iteracao
	do {
		// formula original
		// y[k+1] = (1 - (1 - y[k]^4)^0.25) / (1 + (1 - y[k]^4)^0.25)
		// formula adaptada, mas eh o mesmo resultado
		// y[k+1] = (1 - raiz(raiz(1 - y[k]^4))) / (1 + raiz(raiz(1 - y[k]^4)))

		mpf_pow_ui(t4, y, 4);	// t4 = y[k]^4
		mpf_ui_sub(t4, 1, t4);	// t4 = 1 - y[k]^4
		mpf_sqrt(t4, t4);		// t4 = raiz(1 - y[k]^4)
		mpf_set(t3, t4);		// t3 = raiz(1 - y[k]^4)
		mpf_sqrt(t3, t3);		// t3 = raiz(raiz(1 - y[k]^4))

		mpf_ui_sub(t1, 1, t3);	// t1 = 1 - raiz(raiz(1 - y[k]^4))
		mpf_add_ui(t2, t3, 1);	// t2 = 1 + raiz(raiz(1 - y[k]^4))

		mpf_div(y, t1, t2);		// y[k+1] = (t1) / (t2)


		// formula original
		//a[k+1] = a[k] * (1 + y[k+1])^4 - 2^(2*k + 3) * y[k+1] * (1 + y[k+1] + y[k+1]^2)
		// formula adaptada, mas eh o mesmo resultado
		//a[k+1] = a[k] * (1 + y[k+1])^4 - 2^(2*k + 3) * (y[k+1] + y[k+1]^2 + y[k+1]^3)

		mpf_pow_ui(t4, y, 3);						// t4= y[k+1]^3
		mpf_pow_ui(t3, y, 2);						// t3= y[k+1]^2
		mpf_add(t3, t3, t4);						// t3 = y[k+1]^2 + y[k+1]^3
		mpf_add(t3, y, t3);							// t3 = y[k+1] + y[k+1]^2 + y[k+1]^3

		mpf_set_str(t4, "2", 10);					// t4 = 2
		mpf_pow_ui(t4, t4, (long int) (2*k + 3));	// t4 = 2^(2*k+3)
		mpf_mul(t2, t4, t3);						// t2 = 2^(2*k+3) * (y[k+1] + y[k+1]^2 + y[k+1]^3)

		mpf_add_ui(t4, y, 1);						// t4 = 1 + y[k+1]
		mpf_pow_ui(t4, t4, 4);						// t4 = (1 + y[k+1])^4
		mpf_mul(t1, a, t4);							// t1 = a[k] * (1 + y[k+1])^4

		mpf_sub(a, t1, t2);							// a[k+1] = t1 - t2


		// guarda os valores para verificar se passou de 10 milhoes de casas decimais
		mpf_set(pi_anterior, pi_atual);		// pi_anterior = pi_atual (passa o pi antigo para pi_anterior)
		mpf_ui_div(pi_atual, 1, a);			// pi = 1 / a (atualiza o pi)


		// imprime a iteracao k
		gmp_printf("# k:%d  pi:%.50Ff\n", k, pi_atual);
		

		// proxima iteracao
		k++;

	// comparacao para ver se conseguiu os 10 milhoes de casas exatos
	} while (mpf_cmp(pi_anterior, pi_atual) != 0);

	
	// imprime a quantidade de iteracoes feitas para achar os 10 milhoes de casas decimais do pi
	printf("> Numero total de iteracoes: %d\n", k);


	// limpa a memoria
	mpf_clears(a, y, t1, t2, t3, t4, pi_anterior, pi_atual, NULL);

}



// threads

void *a_initial(void *arg) {
	// pa[0] = 6 - 4 * raiz(2)
	mpf_sqrt_ui(pa, 2);		// pa[0] = raiz(2)
	mpf_mul_ui(pa, pa, 4);	// pa[0] = raiz(2) * 4
	mpf_ui_sub(pa, 6, pa);	// pa[0] = 6 - pa	

	// barreira 1, so libera depois que a outra barreira for executada
	pthread_barrier_wait(&bar1);
	return (NULL);
}

void *y_initial(void *arg) {
	// py[0] = raiz(2) - 1
	mpf_sqrt_ui(py, 2);		// py[0] = raiz(2)
	mpf_sub_ui(py, py, 1);	// py[0] = py[0] - 1

	// barreira 1, so libera depois que a outra barreira for executada
	pthread_barrier_wait(&bar1);
	return (NULL);
}

void *y_numerador(void *arg) {
	mpf_ui_sub(pt1, 1, pt3);	// pt1 = 1 - raiz(raiz(1 - y[k]^4))

	// barreira 2, so libera depois que a outra barreira for executada
	pthread_barrier_wait(&bar2);
	return (NULL);
}

void *y_denominador(void *arg) {
	mpf_add_ui(pt2, pt3, 1);	// pt2 = 1 + raiz(raiz(1 - y[k]^4))

	// barreira 2, so libera depois que a outra barreira for executada
	pthread_barrier_wait(&bar2);
	return (NULL);
}

void *a_division1(void *arg) {
	mpf_pow_ui(pt6, py, 3);		// pt6= py[k+1]^3
	mpf_pow_ui(pt3, py, 2);		// pt3= py[k+1]^2
	mpf_add(pt3, pt3, pt6);		// pt3 = py[k+1]^2 + py[k+1]^3
	mpf_add(pt3, py, pt3);		// pt3 = py[k+1] + py[k+1]^2 + py[k+1]^3

	// barreira 3, so libera depois que a outra barreira for executada
	pthread_barrier_wait(&bar3);
	return (NULL);
}

void *a_division2(void *arg) {
	// pega o valor k passado por argumento
	data *d;
	d = (data *) arg;

	mpf_set_str(pt5, "2", 10);										// t5 = 2
	mpf_pow_ui(pt5, pt5, (long int) (2*(d->number_iterate) + 3));	// t5 = 2^(2*k+3)

	// barreira 3, so libera depois que a outra barreira for executada
	pthread_barrier_wait(&bar3);
	return (NULL);
}

void *a_division3(void *arg) {
	mpf_add_ui(pt4, py, 1);		// pt4 = 1 + py[k+1]
	mpf_pow_ui(pt4, pt4, 4);	// pt4 = (1 + py[k+1])^4
	mpf_mul(pt1, pa, pt4);		// pt1 = pa[k] * (1 + py[k+1])^4

	// barreira 4, so libera depois que a outra barreira for executada
	pthread_barrier_wait(&bar4);
	return (NULL);
}



void calculationPiParBorwein() {

	// variaveis das threads
	pthread_t td1, td2, td3, td4, td5, td6, td7;

	// valores para ser passado por argumentos pras threads
	data d1;

	// valor da iteracao atual
	int k = 0;


	// inicializacao de a, y e temporarios pelo gmp
	mpf_init2(pa, 10000000);
	mpf_init2(py, 10000000);
	mpf_init2(pt1, 10000000);
	mpf_init2(pt2, 10000000);
	mpf_init2(pt3, 10000000);
	mpf_init2(pt4, 10000000);
	mpf_init2(pt5, 10000000);
	mpf_init2(pt6, 10000000);
	mpf_init2(ppi_anterior, 10000000);
	mpf_init2(ppi_atual, 10000000);


	// inicializacao das barreiras de controle do uso das threads
	pthread_barrier_init(&bar1, NULL, 2);
	pthread_barrier_init(&bar2, NULL, 2);
	pthread_barrier_init(&bar3, NULL, 2);
	pthread_barrier_init(&bar4, NULL, 2);


	// valores iniciais
	// cria e executa as threads t1 e t2
	pthread_create(&td1, NULL, (void *) &a_initial, NULL);
	pthread_create(&td2, NULL, (void *) &y_initial, NULL);

	// junta as threads t1 e t2
	pthread_join(td1, NULL);
	pthread_join(td2, NULL);


	// iteracao
	do {
		// formula original
		// y[k+1] = (1 - (1 - y[k]^4)^0.25) / (1 + (1 - y[k]^4)^0.25)
		// formula adaptada, mas eh o mesmo resultado
		// y[k+1] = (1 - raiz(raiz(1 - y[k]^4))) / (1 + raiz(raiz(1 - y[k]^4)))

		mpf_pow_ui(pt4, py, 4);		// t4 = y[k]^4
		mpf_ui_sub(pt4, 1, pt4);	// t4 = 1 - y[k]^4
		mpf_sqrt(pt4, pt4);			// t4 = raiz(1 - y[k]^4)
		mpf_set(pt3, pt4);			// t3 = raiz(1 - y[k]^4)
		mpf_sqrt(pt3, pt3);			// t3 = raiz(raiz(1 - y[k]^4))

		// cria as threads t3 e t4
		pthread_create(&td3, NULL, (void *) &y_numerador, NULL);
		pthread_create(&td4, NULL, (void *) &y_denominador, NULL);

		// junta as threads t3 e t4
		pthread_join(td3, NULL);
		pthread_join(td4, NULL);

		mpf_div(py, pt1, pt2);		// y[k+1] = (t1) / (t2)


		// formula original
		//a[k+1] = a[k] * (1 + y[k+1])^4 - 2^(2*k + 3) * y[k+1] * (1 + y[k+1] + y[k+1]^2)
		// formula adaptada, mas eh o mesmo resultado
		//a[k+1] = a[k] * (1 + y[k+1])^4 - 2^(2*k + 3) * (y[k+1] + y[k+1]^2 + y[k+1]^3)

		// passa o numero de iteracoes para uma thread
		d1.number_iterate = k;

		//cria as threads t5, t6 e t7
		pthread_create(&td5, NULL, (void *) &a_division1, NULL);
		pthread_create(&td6, NULL, (void *) &a_division2,  (void *) &d1);
		pthread_create(&td7, NULL, (void *) &a_division3, NULL);

		// junta as threads t5 e t6
		pthread_join(td5, NULL);
		pthread_join(td6, NULL);

		mpf_mul(pt2, pt5, pt3);		// t2 = 2^(2*k+3) * (y[k+1] + y[k+1]^2 + y[k+1]^3)

		// barreira para finalizar a td7 ou o calculo acima		
		pthread_barrier_wait(&bar4);

		// junta a thread t7 com o calculo acima
		pthread_join(td7, NULL);

		mpf_sub(pa, pt1, pt2);		// a[k+1] = t1 - t2


		// guarda os valores para verificar se passou de 10 milhoes de casas decimais
		mpf_set(ppi_anterior, ppi_atual);	// pi_anterior = pi_atual (passa o pi antigo para pi_anterior)
		mpf_ui_div(ppi_atual, 1, pa);		// pi = 1 / a (atualiza o pi)


		// imprime a iteracao k
		gmp_printf("# k:%d  pi:%.50Ff\n", k, ppi_atual);
		

		// proxima iteracao
		k++;


	// comparacao para ver se conseguiu os 10 milhoes de casas exatos
	} while (mpf_cmp(ppi_anterior, ppi_atual) != 0);

	
	// imprime a quantidade de iteracoes feitas para achar os 10 milhoes de casas decimais do pi
	printf("> Numero total de iteracoes: %d\n", k);

	// limpa a memoria
	mpf_clears(pa, py, pt1, pt2, pt3, pt4, pt5, pt6, ppi_anterior, ppi_atual, NULL);
}
